/*
 * fft.c
 *
 *  Created on: 2009-7-24
 *      Author: wyb
 */

#include "fft.h"

void output_complex(struct fft_complex *value, char *str) {
	printf("--------------------\n");
	printf("%s\n", str);
	printf("--------------------\n");
	if (value == NULL) {
		printf("NULL");
		return;
	}
	int i;
	for (i = 0; i < value->N; ++i) {
		printf("%.4f %.4f \n", value->data[i][0], value->data[i][1]);
	}
	printf("--------------------\n");
}

void file_complex(struct fft_complex *value, char *str) {
	FILE *fp;
	if ((fp = fopen(str, "w+")) == NULL) {
		printf("can't open %s", str);
	}

	int i;
	for (i = 0; i < value->N - 1; i++) {
		fprintf(fp, "%f,", value->data[i][0]);
	}
	fprintf(fp, "%f", value->data[value->N - 1][0]);
	fclose(fp);
}

void output_cos(struct fft_cosine *value, char *str) {
	printf("--------------------\n");
	printf("%s\n", str);
	printf("--------------------\n");
	if (value == NULL) {
		printf("NULL");
		return;
	}
	int i;
	for (i = 0; i < value->N; ++i) {
		printf("%.4f %.4f %.4f \n", value->A[i][0], value->A[i][1],
				value->A[i][2]);
	}
	printf("--------------------\n");
}

void file_cos_A(struct fft_cosine *value, char *str) {
	FILE *fp;
	if ((fp = fopen(str, "w+")) == NULL) {
		printf("can't open %s", str);
	}

	int i;
	for (i = 0; i < value->N - 1; i++) {
		fprintf(fp, "%f,", value->A[i][0]);
	}
	fprintf(fp, "%f", value->A[value->N - 1][0]);
	fclose(fp);
}

struct fft_complex *fft(struct fft_complex *value) {
	struct fft_complex *result;
	fftw_plan plan;
	result = malloc(sizeof(struct fft_complex));
	result->N = value->N;
	result->data = (fftw_complex*) fftw_malloc(value->N * sizeof(fftw_complex));

	plan = fftw_plan_dft_1d(value->N, value->data, result->data, FFTW_FORWARD,
			FFTW_ESTIMATE);

	fftw_execute(plan);

	fftw_destroy_plan(plan);

	return result;
}

struct fft_complex *ifft(struct fft_complex *value) {
	struct fft_complex *result;
	fftw_plan plan;
	result = malloc(sizeof(struct fft_complex));
	result->N = value->N;
	result->data = (fftw_complex*) fftw_malloc(value->N * sizeof(fftw_complex));

	plan = fftw_plan_dft_1d(value->N, value->data, result->data, FFTW_BACKWARD,
			FFTW_ESTIMATE);

	fftw_execute(plan);

	fftw_destroy_plan(plan);

	int i;
	for (i = 0; i < result->N; i++) {
		result->data[i][0] = result->data[i][0] / result->N;
	}

	return result;
}

struct fft_complex *ifft_cosine(struct fft_cosine *value, int n) {
	struct fft_complex *result;
	result = malloc(sizeof(struct fft_complex));
	result->N = n;
	result->data = (fftw_complex*) fftw_malloc(n * sizeof(fftw_complex));

	int i, j;
	for (i = 0; i < n; i++) {
		double A = 0;
		for (j = 0; j < value->N; j++) {
			A += value->A[j][0] * cos(2 * PI * value->A[j][2] * i
					+ value->A[j][1]);
		}
		result->data[i][0] = A;
		result->data[i][1] = 0.0;
	}

	return result;
}

double compute_conine(struct fft_cosine *value, int N, double period, double x) {
	double x_temp = x - ((int) (x / period)) * period;
	int i, j;
	double result = 0.0;
	for (i = 0; (x_temp + i * period) <= (N - 1); i++) {
		double temp = 0.0;
		double y = x_temp + i * period;
		for (j = 0; j < value->N; j++) {
			temp += value->A[j][0] * cos(2 * PI * value->A[j][2] * y
					+ value->A[j][1]);
		}
		result += temp;
	}
	return (result / i);
}

struct fft_cosine *complex_to_cosine(struct fft_complex *value) {
	struct fft_cosine *result;
	result = malloc(sizeof(struct fft_cosine));
	result->N = value->N;
	result->A = malloc(value->N * sizeof(double[3]));

	int i;
	for (i = 0; i < value->N; i++) {
		double a = value->data[i][0];
		double b = value->data[i][1];
		result->A[i][0] = sqrt(a * a + b * b) / value->N;
		result->A[i][1] = atan2(b, a);
		result->A[i][2] = (double) i / value->N;
	}

	return result;
}

struct fft_complex *real_to_complex(struct double_array *value) {

	struct fft_complex *result;
	result = malloc(sizeof(struct fft_complex));
	result->N = value->N;
	result->data
			= (fftw_complex*) fftw_malloc(result->N * sizeof(fftw_complex));

	int i;
	for (i = 0; i < result->N; i++) {
		result->data[i][0] = value->data[i];
		result->data[i][1] = 0.0;
	}

	return result;
}

struct double_array *complex_to_real(struct fft_complex *value) {
	struct double_array *result;
	result = malloc(sizeof(struct double_array));
	result->N = value->N;
	result->data = malloc(value->N * sizeof(double));

	int i;
	for (i = 0; i < value->N; i++) {
		result->data[i] = value->data[i][0];
	}

	return result;
}
